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Abstract Ecological modelling of increasingly more 
complex microbial populations is necessary to reflect 
the highly functional and diverse behaviour inherent 
to many systems found in reality. Anaerobic digestion 
is one such process that has benefitted from the ap¬ 
plication of mathematical analysis not only for char¬ 
acterising the biological dynamics, but also to inves¬ 
tigate emergent behaviour not apparent by simulation 
alone. Nevertheless, the standard modelling approach 
has been to describe biological systems using sets of 
differential equations whose kinetics are generally de¬ 
scribed by some empirically derived function of growth. 
The drawbacks of this are two-fold; the growth func¬ 
tions are derived from empirical studies that may not 
be representative of the system to be modelled and 
whose parameters may not have a mechanistic mean¬ 
ing, and mathematical analysis is restricted by a con¬ 
formity to an assumption of the dynamics. Here, we 
attempt to address these challenges by investigating a 
generalised form of a three-tier chlorophenol mineral¬ 
ising food-web previously only analysed numerically. 
We examine the existence and stability of the identi¬ 
fied steady-states and find that, without a decay term, 
the system may be characterised analytically. However, 
it is necessary to perform numerical analysis for the 
case when maintenance is included, but in both cases 
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we verify the discovery of two important phenomena; 
i) the washout steady-state is always stable, and ii) the 
two other steady-states can be unstable according to 
the initial conditions and operating parameters. 

Keywords Microbial modelling ■ Dynamical systems • 
Stability theory • Anaerobic digestion 

1 Introduction 

The mathematical modelling of engineered biological 
systems has entered a new era in recent years with the 
expansion and standardisation of existing models aimed 
at collating disparate components of these processes 
and provide scientists, engineers and practitioners with 
the tools to better predict, control and optimise them. 
These forms of mechanistic models emerged initially 
with the Activated Sludge Models m for wastewater 
treatment processes, followed by the Anaerobic Diges¬ 
tion Model No. 1 (ADM1) |3] a few years later. The 
development of ADM1 was enabled largely due to the 
possibilities for better identification and characterisa¬ 
tion of functional groups responsible for the discrete 
degradation steps operating in series within anaerobic 
digesters. It describes a set of fairly complex stoichio¬ 
metric and kinetic functions representing the standard 
anaerobic process, remaining the scientific benchmark 
to the present day, despite an general understanding of 
its limitations in describing all necessary biochemical 
transformations. Indeed, there has been a growing argu¬ 
ment that the model should take advantage of improved 
empirical understanding and extension of biochemical 
processes included in its structure, to acquire a better 
trade-off between model realism and complexity [31 . 

It has previously been shown that the simplification 
or reduction in model complexity can preserve biolog- 
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ical meaning whilst reducing the computational effort 
required to find mathematical solutions of the model 
equations mm- Whilst simpler models are approx¬ 
imations of real systems, it can be beneficial to con¬ 
sider a reduced model to better understand biologi¬ 
cal phenomena of sub-processes without the need to 
consider extraneous system parameters and variables, 
which tend to make mathematical analysis intractable 
and cumbersome. Nevertheless, even with gross simpli¬ 
fication of a biological system based on a set of ordinary 
differential equations (ODEs) of relatively low dimen¬ 
sionality, analytical techniques are unable to provide 
general solutions for the system and numerical meth¬ 
ods must suffice. 

As an example of this, for anaerobic digestion, a 
previous study investigated the effect of maintenance 
on the stability of a two-tiered ‘food-chain’ compris¬ 
ing two species and two substrates J5]. Although the 
authors were not able to determine the general con¬ 
ditions under which this four dimensional syntrophic 
consortium was stable, further work has shown that a 
model with generality can be used to answer the ques¬ 
tion posed, determining that the two-tiered food-chain 
is always stable when maintenance is included [5]. 

In a more recent example, the model described by [5] 
was extended by the addition of a third organism and 
substrate to create a three-tiered ‘food-web’ [lD]. In 
this model, the stability of some steady-states could be 
determined analytically, but due to the complexity of 
the Jacobian matrix for certain steady-states, local so¬ 
lutions were necessary using numerical analysis, when 
considering the full system behaviour. Although the re¬ 
sults were important in revealing emergent properties 
of this extended model, the motivation of this work is 
to determine whether the approach carried out in [5], 
can be applied to the three-tiered model from m , to 
provide some general properties of that system. 

The paper is organised as follows. In Section [2] we 
present a description of the model to be investigated, 
before providing an alternative reduction of its struc¬ 
ture than that given by nni, in Section [3j With Sec- 
tion[4]we demonstrate the existence of the three steady- 
states and define four interesting cases for specific pa¬ 
rameter values that are investigated using the analytical 
solutions, whilst also indicating the regions of existence 
of the steady-states for the operating parameter val¬ 
ues (dilution rate and substrate input concentration). 
In Section [5] we perform local stability analysis of the 
steady-states without maintenance, and in Section [6j 
perform a comprehensive numerical stability analysis 
of the four cases for both the model with and without a 
decay constant. We show that our approach leads to the 
discovery of five operating regions, in which one leads 


to the possibility of instability of the positive steady 
state, where all three organisms exist, a fact that has 
not be reported by ina- Indeed, we prove that a stable 
limit-cycle can occur in this region. Finally, in Section]?} 
we make comment on the role of the kinetic parameters 
used in the four example cases, in maintaining stability, 
which points to the importance of the relative aptitude 
of the two hydrogen consumers in sustaining a viable 
chlorophenol mineralising community. In the Appendix 
we describe the numerical method used in Section[6]and 
we give the proofs of the results. 


2 The model 

The model developed in nni has six components, three 
substrate (chlorophenol, phenol and hydrogen) and 
three biomass (chlorophenol, phenol and hydrogen de¬ 
graders) variables. The substrate and biomass concen¬ 
trations evolve according to the six-dimensional dynam¬ 
ical of ODEs 


dA, 


ch 


d t 

dXph 
d t 

d.V H2 
d t 

dSch 

di 

dS ph 


= ~DX C h + Ytfifo (Sell, Sh 2 ) A' ch - fcdec.chWh (1) 

= —DXp h + Yph/l (Sph, Sh 2 ) Xph — fcdec,phW>h (2) 

= —DX h 2 +}h 2 /2 (Sh 2 )Xh 2 — fcdec,H 2 ^H 2 (3) 

= D (S ch , in - S ch ) - fo (S ch , Sh 2 ) x ch (4) 

224 


dt = D (S ph , in - s ph ) + — (1 - y ch ) fo (S ch , Sh 2 ) A' ch 
-/i (S ph ,S H2 )A' ph (5) 

= (Sh 2 ,in — Sh 2 ) + (1 — Wh) fl (Sph, Sh 2 ) A'ph 


16 

208 


fo (S c h, Sh 2 ) A ch - A (Sh 2 ) A"h 2 


( 6 ) 


where S c h and X c h are the chlorophenol substrate and 
biomass concentrations, S p h and X p h those for phe¬ 
nol and Sua and A'h 2 those for hydrogen; Y c h, F p h 
and Yh 2 are the yield coefficients, 224/208 (1 — F c h) 
represents the part of chlorophenol degraded to phe¬ 
nol, and 32/224(1 — Yph) represents the part of phe¬ 
nol that is transformed to hydrogen. Growth functions 
take Monod form with hydrogen inhibition acting on 
the phenol degrader and represented in fi (see Eq. [7| 
as a product inhibition term. 


f ( C Q \ _ fcm.chSch S/J 2 

Jo Wch, oh 2 J - Ks,oh+S ch K 3lM2 .f+Sm 2 

fi {S ph ,Sii 2 ) = i-Ho (7) 

+ K i,H 2 

r / C \ km,H 2 S H 2 

J2 (*H 2 j — Ks,„ 2 *S„ 2 

Here, apart from the four operating (or control) pa¬ 
rameters, which are the inflowing concentrations .S/h.in, 
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5ph,in ; *5h2 ,m toid the dilution rate D , that can vary, all 
others have biological meaning and are fixed depending 
on the organisms and substrate considered. We use the 
following simplified notations in (Eqs. [lj6]) 


Xq = X ch , X, = X ph , = Xh 2 
Sq O c h - *01 — ■ ^2 b'l [,, 

cin _ c cin _ O Qin _ q 

*^0 — *^ch,in; *^1 — ^ph.ini ^2 — *oH 2 ,ir] 

Y 0 = Fch, 

224 


ii = Iph, 


Y2 = Wi 2 

y 4 = — (1 - y ph ), y 5 = — 

224 ' p 5 208 


Ys 208 ^ Fch ^ 1 
ItO — ^-dec.ch, til — ^-dec.ph, M 2 — ^dec,H 2 

With these notations Eqs.[l](6]can be written as follows 


We obtain the following system 


dxo 

dt 

= -Dxo 

+ MO 

(so, s 2 ) rro 

— CLOXO 

(14) 

cirri 

dt 

= —Dx 1 

+ m 

(si,s 2 ) XI 

— a\x\ 

(15) 

dx 2 

dt 

= —Dx 2 

+ M 2 

(S 2 ) x 2 ~ a 2 X 2 

(16) 

dso 

dt 

= D (sq 11 

- so 

O 

O 

3 . 

1 

82)3:0 

(17) 

ds\ 

dt 

= —Ds 1 

+ Mo 

(so,s 2 )x 0 

- Ml (81,82)3:1 

(18) 

ds 2 

dt 

= —Ds 2 

+ Mi 

(si,S 2 ) xi 

- Uipo (so, 8 2 ) 3:0 

- M2 (82)3:2 





(19) 


where the inflowing concentration is 


dX 0 

dt 

dXi 

dt 

dX 2 

dt 

dSp 

dt 

dSi 

dt 

dS 2 

dt 


= -DX 0 + Y 0 f 0 (So,S 2 )Xo-a 0 X 0 ( 8 ) 

= -DX 1 + Y 1 f 1 (S 1 ,S 2 )X 1 -a 1 X 1 (9) 

= -DX 2 + Y 2 f 2 (S 2 ) X 2 - a 2 X 2 (10) 

= D(S i o n -S o )-fo(So,S 2 )X 0 (11) 

= D (S[ n - Si) + Y 3 f 0 (So, S 2 ) X 0 - fi (Si,S 2 ) Xi (12) 

= D (S? - S 2 ) + > 4/1 (Si, S 2 ) Xi - y 5 / 0 (So, S 2 ) X 0 
-/ 2 (S 2 )X 2 (13) 


In (TO], this model is reduced to a dimensionless form 
that significantly reduces the number of parameters de¬ 
scribing the dynamics. In this paper we do not assume 
that the growth functions /o, fi and f 2 have the specific 
analytical expression (Eq. [?| . We will only assume that 
the growth functions satisfy properties that are listed 
in Appendix [Cj Therefore, we cannot benefit from the 
dimensionless rescaling used by [10], because this rescal¬ 
ing uses some kinetics parameters of the specific growth 
functions (Eq. [t]) , while we work with general unspeci¬ 
fied, growth functions. In Section [3] we consider another 
rescaling that does not use the kinetics parameters. Fur¬ 
thermore, we restrict our analysis to the case where we 
only have one substrate addition to the system, such 
that: Sifi > 0, S\ n = 0, and S l 2 n = 0. 


3 Model reduction 


To ease the mathematical analysis, we can rescale the 
system (Eqs. 8|TTT|) using the following change of vari¬ 
ables adapted from [9] : 


I 3 I 4 „ 

Xo — v A 0 , 

to 

SO = Y3Y4S0, 


Xl = wA 1 , 

fl 


12 = 


Si — Y 4 S\, s 2 — S 2 


sir = >3 w, 

the growth functions are 

A®o(soj s 2 ) = k 0 /o ^ Y 3 Y 4 j 
Mi(si,s 2 ) = W/i (ft>s 2 ) 
A* 2 (s 2 ) = I 2 / 2 (s 2 ) 
and 

= JL = 1 

w y 3 i 4 2(1 - y 0 )(i - y x ) 


( 20 ) 


( 21 ) 


( 22 ) 


The benefit of our rescaling is that it permits to fix 
in Eqs. [T4]|T9| all yield coefficients to one except that 
denoted by u and defined by (Eq. 22), and to discuss 
the existence and stability with respect to this sole pa¬ 
rameter. 

Using Eq. 21 and the growth functions (Eq. [Tf, 
we obtain the model (Eqs. p~T|To[ ) with the following 
Monod-type growth functions 


Mo (so, s 2 ) 
Mi (si,s 2 ) 


mpsp s 2 
Kp+sp Lp+S 2 


m\S \1 


M 2 (s 2 ) 


m 2 S 2 

K 2 +S 2 


where 


(23) 


UlO — ^O^m.ch, Kq — Ao — IlS,H 2) c 

I^l^mjphj A| Ul-I’Cs.ph, X-i (24) 

m 2 = y 2 fc m ,H 2 , I\ 2 = K Ts,h 2 

For the numerical simulations we will use the nominal 
values in Table [l] given in m ■ 
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Parameters 

Nominal values 

Units 

^m,ch 

29 

kgCODs/kgCODx/d 

Ks, ch 

0.053 

kgCOD/m 3 

Uch 

0.019 

kgCODx/kgCOD s 

^m,ph 

26 

kgCODs/kgCODx/d 

K S,ph 

0.302 

kgCOD/m 3 

>ph 

0.04 

kgCODx/kgCODs 

^m, H2 

35 

kgCODs /kgCODx/d 

a 's,h 2 

2.5x 10 -5 

kgCOD/m 3 

^S,H 2 ,c 

l.Ox 10 -6 

kgCOD/m 3 

Vh 2 

0.06 

kgCODx/kgCODs 

^dec,i 

0.02 

d” 1 

Ki, h 2 

3.5x 10 -6 

kgCOD/m 3 


Table 1 Nominal parameter values. 


4 Existence of steady-states 

A steady-state of Eqs. PP 1 is obtained by setting the 
right-hand sides equal to zero: 


[mo (s 0 ,s 2 ) - D - a 0 \x 0 = 0 (25) 

[mi (si, s 2 ) - D - ai] xi = 0 (26) 

[M2 (s 2 ) - D - a 2 ] x-2 = 0 (27) 

D (sg 11 — s 0 ) — fio (so,S 2 )xo = 0 (28) 

—Ds i + /to (so,S 2 )xo — Mi (si, S 2 ) xi = 0 (29) 

-Ds 2 + Mi (si,s 2 )x 1 — oj/jo ( s 0 ,s 2 )x 0 - M 2 (s 2 )x 2 = 0 (30) 

A steady-state exists (or is said to be ‘meaningful’) if, 
and only if, all its components are non-negative. 

Lemma 1 The only steady-state of Eqs. ~L$T9 for 
which xq = 0 or x\ = 0, is the steady-state 

SSI = (x 0 = 0, xi = 0, x 2 = 0, s 0 = s™, si = 0,s 2 = 0) 


where all species are washed out. This steady-state al¬ 
ways exists. It is always stable. 


From the previous Lemma we deduce that besides 
the steady-state SSI, the system can have at most two 
other steady-states. 

SS2: xo > 0, xi > 0 and x 2 = 0, where species x 2 is 
washed out while species xo and and xi exist. 

SS3: Xo > 0, Xi > 0, and x 2 > 0, where all popula¬ 
tions are maintained. 

In the following we describe the steady-states SS2 
and SS3 of Eqs. TlfTT)] with the Monod-type growth 
functions (Eq. [23 1. The general case with unspecified 
growth function is provided in Appendix [Cj with proofs 
given in Appendix [Dj We use the following notations: 

Let s 2 be fixed, we define the function Mo(y, s 2 ) as 


follows : for all y £ 


0,/xo(+oo,s 2 ) 


moS2 

Lq + S2 


we let 


M 0 (y,s 2 ) 


K 0 y 


m 0 s 2 _ 
L 0 +S 2 " 


Notice that y i—>■ Mo(y, s 2 ) is the inverse function of the 
function So H>• yo(so,s 2 ), that is to say, for all sq > 0, 
s 2 > 0 and y £ [0, y 0 (+oo, s 2 )) 

s 0 = M 0 (y, s 2 ) y = Mo(so,s 2 ) (31) 


Let s 2 be hxed, we define the function Mi(y, s 2 ) as 
follows : for all y £ 0, /zi(+oo, s 2 ) = i +2)k) i we let 


M 1 (y,s 2 ) = 


Kiy 


m i 

1 +s 2 /Ki 


Notice that y i—>■ Mi(y, s 2 ) is the inverse function of the 
function si i—>• Mi(si,s 2 ), that is to say, for all Si > 0, 
s 2 > 0 and y £ [0,/xi(+oo, s 2 )) 


si = Mi(y, s 2 ) y = Mi(si, s 2 ) 


(32) 


We define the function M 2 (s 2 ) as follows : for all 
y £ [0, /i 2 (Too) = m 2 ), we let 


M 2 (y) = 


I< 2 V 
m 2 — y 


Notice that y i—>• M 2 (y) is the inverse function of the 
function s 2 i— > y 2 {s 2 ), that is to say, for all s 2 > 0 and 
V e [0,/Li 2 (+oo)) 


s 2 = M 2 (s 2 ) 


y = M 2 (s 2 ) 


(33) 


Using the functions Mo, Mi and M 2 we dehne the 
following function: Let u> < 1. Let 

\ rn i \ i + ai, s 2 ) + s 2 

ip{s 2 ) = M 0 (D+ ao,s 2 ) -1 -t- ( 34 ) 

l — UJ 

Notice that if is defined if, and only if, 

D + a 0 < jj, o(+oo, s 2 ) and D + a\ < Mi(+oo, s 2 ) 


which is equivalent to 

s°(D) < s 2 < s\(D) 


where 


1(D) = 


Lq (D + ao) 


m o — D — ao 
are the solutions of equations 


1(D) 


Kj(mi - D - ai) 
D + ai 


Mo(+oo, s 2 ) = D + ao, pi(+oo, s 2 ) = D + a\ (35) 

respectively, see Fig. [I] (a). Straightforward calculations 
show that 

_ T K 1 (K i +s 2 ) , 0 

,f i —, \ Dq(D + ao) Lq + s 2 s\(D)—s 2 ' 2 

if(s 2 ,D) = --- o7Ty\ + - 

/no — U — ao s 2 — S 2 (T>) 1 — co 

Therefore, if(s 2 ) > 0 for s° < s 2 < S 2 (see Fig. [l](b)), 


lim if(s 2 ) = lim if(s 2 ) = +00 

S2—S2—>-<82 
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Fig. 2 Graphs of an< i s h(D) (' n black) and M 2 (D) (in red) and graphical depiction of I\ = [0, Di), where D\ is the 

solution of s 2 (D) = s\(D), and I 2 ■ (a): I 2 = [0, Z? 2 ) where D 2 is the solution of M 2 (D) = s\(D). (b) : I 2 = [ 0 , 1 ) 2 ) where D 2 
is the solution of M 2 (D) = s 2 (D). (c): I 2 is empty, (d) : I 2 = (D 2 min, ftmoi) where D 2 min and D 2 max are the solutions of 
M 2 (D) = s 2 (D) and M 2 (D ) = s}(D), respectively. 


and 

d 2 ij) = 2K 0 (D + a 0 ) L 0 + s°{D) _ 2K 1 {K i + 4(D)) 
ds\ ~ mo- D- a 0 ( S2 - s°(D)) 3 (1 - u) (4(D) - s 2 ) 3 

Hence, > 0 for all s 2 G (s^D), s 2 (D)), so that the 
function S 2 1 —> ij)(s 2 ,D) is convex and, thus, it has a 
unique minimum 132 (D), see Fig. 0(b)- 
Let w < 1. We define the function 

F 1 (D)= in| ip(s 2 ) = ip (s 2 ) (36) 

s 2 e(s§,s^) 

as shown in Fig. 0 (b). The minimum s 2 (D) is a so¬ 
lution of an algebraic equation of degree 4 in s 2 . Al¬ 
though mathematical software, such as Maple , cannot 
give its solutions explicitly with respect to the param¬ 
eters, s 2 (D) could be obtained analytically since alge¬ 
braic equations of degree 4 can theoretically be solved 
by quadratures. We do not try to obtain such an explicit 
formula. However, if the biological parameters are fixed, 
the function s 2 (D) and, hence, F\(D) = tp(s 2 (D), D), 
can be obtained numerically. 

The function F\(D) is defined as long as s 2 (D) < 
4(D). Assuming that 4(0) < 4(0), Fi(D) is defined 
for 0 < D < Di, where D\ is the positive solution of 
s° 2 (D)=s\(D) (see Fig. § . Therefore, D\ is a solution 
of the the second order algebraic equation = 


Ki( ' ni D + a 1 ° 1 ' 1 . We denote by 

h = {D: 4(D) < 4(D)} (37) 

the set on which F\(D) is defined. 

Let w < 1. We define the functions 


F 2 (D) = ^ (M 2 (D + a 2 )) 
F 3 (D) = (M 2 (D + a 2 )) 

CIS 2 


(38) 

(39) 


Since M 2 and xjj are given explicitly by Eq. [33] and 


Eq. 34 respectively, the functions F 2 (D) and F 3 (D) 


are given explicitly with respect to the biological pa¬ 


rameters in Eq. 23 The functions F 2 (D) and F 3 (D) 
are defined for D such that s 2 (D) < M 2 (D) < s\(D ), 
that is to say, for D such that 

L 0 (D + a 0 ) K 2 (D + a 2 ) Ki(m\ — D — a{) 


mo — D — ao 
We denote by 


< 


m 2 — D — a 2 


< 


D + a± 


I 2 = {De /, : 4(D) < M 2 (D) < 4(D)} 


(40) 


the subset of I\ on which F 2 (D) and F 3 (D) are de¬ 
fined. For all for D G I 2 , F\(D) < F 2 (D). The equal¬ 
ity F\(D) = F 2 (D) holds if, and only if, M 2 (D + 
a 2 ) = s 2 (D) that is, ^ (M 2 (D + a 2 )) = 0. Therefore, 
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Fi(D) = F 2 (D) holds if, and only if, F 3 (D) = 0. We gives a steady-state SS2 = (xo,x 4 ,x 2 = 0, So, si, S 2 ) 
define where 


I 3 = {D € I 2 : F 3 (D) < 0} 

Since D i-> s 2 (D) is increasing and D 1 —>■ s 2 (D) is 
decreasing, and assuming s 2 (0) < s 2 (0), the domain 
of definition I\ of F\(D) is an interval I\ = [0, D 4 ), 
where D\ is the solution of s 2 (D) = s\(D), see Fig. [ 2 ] A 
necessary condition of existence of SS2 is 0 < D < D\. 

For the domain of definition J 2 of F 2 (D ), several 
cases can be distinguished. I 2 is an interval I 2 = [0, D 2 ), 
where D 2 is the solution of M 2 (D) = s\(D), see Fig. 
H a )> or the solution of equation M 2 (D) = So(D), see 
Fig.§b). h is empty, see Fig. Ha)- I 2 is an interval 
I 2 = ( D 2 m i ni D 2max ) where t) 2 rll j /n and D 2max are the 
solutions of M 2 {D) = s 2 (D) and M 2 {D) = s\(D) re¬ 
spectively, see Fig. Hd). A necessary condition of exis¬ 
tence of SS3 is D £ I 2 . Cases (a)-(d) are obtained with 
the numerical parameter values listed in Table [2] and [3] 



^S,H 2 ,c 

CLi 

Di 

d 2 

d 3 

(a) 

i.o x io ~ 6 

0.02 

0.432 

0.373 

0.058 



0 

0.452 

0.393 

0.078 

(b) 

4.0 x 10 -6 

0.02 

0.329 

0.236 

h = h 



0 

0.349 

0.256 

A 1 

II 

(c) 

7.0 x 10 -6 

0.02 

0.287 

I 2 ”0 




0 

0.303 

Is = 0 



Table 2 Parameter values for cases (a), (b) and (c) of Fig. 
[ 2 ] Unspecified parameter values are as in Table [I] The table 
gives the values of D\, D 2 and D 3 where I\ = [0, Di), I 2 = 
[ 0 ,D 2 ) and I 3 = [ 0 ,D 3 ) 


di 

Di 


D 2 max 

d 3 

(d) 0.02 

0.238 

0.101 

0.198 

0.161 

0 

0.258 

0.121 

0.218 

0.181 


Table 3 Parameter values for case (d) of Fig. IHI K S , H 2 ,c — 
1.2 x 10 -5 , Ks, h, = 0.5 x 10 -5 and k m: h 2 = 5. Unspeci¬ 
fied parameter values are as in Table ^ The table gives the 
values of D 1 , D 2 min, D 2 max and D 3 where h = [0,Di), 
I 2 = ( 02 mmj 02 ma i) and I 3 = (D 2 min , U 3 ). 


We can state now the necessary and sufficient con¬ 
ditions of existence of SS2 and SS3. 


Lemma 2 If oj > 1 then SS2 does not exist. If to < 1 
then SS2 exists if, and only if, Sq" > F\(D). Therefore, 
a necessary condition for the existence of SS2 is that 
D £ I\, where I\ is defined by Eq. 37. If Sg n > F\(D) 
then each solution s 2 of equation 


V’(s 2 ) = sJ ) n , s 2 e(s^,4) 


(41) 


So — Mo(D + ao, s 2 ), Si — Mi(D + a i, s 2 ) 

a-’o yy I 7 (^o ^o), x 4 yy j ~ (^g So si) 

(42) 


D 


a 0 


D -\- a\ 


Lemma 3 If uj > 1 then SS3 does not exist. If w < 1 
then SS3 exists if, and only if, s™ > F 2 (D). Therefore, 
a necessary condition of existence of SS3 is that D £ I 2 , 
where I 2 is defined by Eq. JO If sf} 1 > F 2 (D) then the 
steady-state SS3 = (xq, #i, x 2 , Sg, si, s 2 ) is given by 


s 0 = M 0 (D + a 0 , M 2 (D + a 2 )) 
si = + ai, M 2 (D + d 2 )) 

s 2 = M 2 (D + a 2 ) (43) 


and 


D 


-Oo n -s 0 ), Xi = 


D 


D + di' 


X ° - 25 + 0 ^°° 

X 2 = —^ ((1 - W)(s(, n - Sg) - Si - S 2 ) 

u + a 2 


(s“ - s 0 - Si) 
(44) 


Remark 1 If s™ > F 4 (D) then Eq. 41 has exactly two 
solutions denoted by s 2 and s 2 and such that, see Fig. 

EH 


s 2 < s 2 < s 2 < s 2 < s 2 

l> _ 


If s“ = F 0 (D) then s 2 < s 2 = s 2 = s 2 < s 2 . 


To these solutions, s 2 and s 2 , correspond two 
steady-states of SS2, which are denoted by SS2 b and 
SS2 I) . These steady-states coalesce when Sg n = F 0 (D). 

Since F\{D) < F 2 (D), the condition s™ > F 2 {D) 
for the existence of the positive steady-state SS3 implies 
that the condition s™ > F 2 (D) for the existence of the 
two steady-states SS2 b and SS2 1 * is satisfied. Therefore, 
if SS3 exists then SS2 b and SS2** exist and are distinct. If 
Sg 11 = F 2 (D) then SS3 coalesces with SS2 b if F 3 (D) < 0, 
and with SS2 # if F 3 (D) > 0, respectively. 


Remark 2 Using Eq. 20 the conditions s™ > F 4 (D) 
and s™ > F 2 (D) of existence of the steady-state SS2 
and SS3 respectively are equivalent to the conditions 


'-’cli.in U 


Fi(D) 

y 3 y 4 


and bVh.in U 


F 2 (D) 

y 3 y 4 


respectively, expressed with respect to the inflowing 
concentration 5 c h,i n - 


Our aim now is to describe the operating diagram : 
The operating diagram shows how the system behaves 
when we vary the two control parameters ^ch.in and D 
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in Eqs. [l][ 6 ] According to Remark [2j the curve F\ of 
equation 

<S'ch I in = yZY^ 1 ^) (45) 

is the border to which SS2 exists, and the curve F 2 of 
equation 

Schdn = YW4 F2{D) (46) 

is the border to which SS3 exists, see Fig. [3] If we want 
to plot the operating diagram we must fix the values 
of the biological parameters. In the remainder of the 
Section we plot the operating diagrams corresponding 
to cases (a)-(d) depicted in Fig. [2j 





Fig. 3 The curves A (black), A (red) and A (green) for case 
(a), (i) : regions of steady-state existence, with maintenance. 
On the right, a magnification for 0 < D < D$ = 0.058 showing 
the region Ji. (ii) : regions of steady-state existence and their 
stability, without maintenance. On the right, a magnification 
for 0 < D < D 3 = 0.078 showing the regions J 4 and J$. 


4.1 Operating diagram: case (a) 

This case corresponds to the parameter values used 
by [10]. We have seen in Table [ 2 ] that the curves A 
and r 2 are defined for D < D\ and D < D 2 , respec¬ 
tively and that they are tangent for D = A 3 , where 
Di = 0.432, D 2 = 0.373 and A 3 = 0.058. Therefore, 
they separate the operating plane {S c ^- ln ,D) into four 
regions, as shown in Fig. [3ji) , labelled J \, J 2 and J 3 
and J 4 . 

The results are summarised in Table [4j which shows 
the existence of the steady-states SSI, SS2 and SS3 in 
the regions of the operating diagram in Fig. [3ji) . 


Region 

Steady states 

Ji 

SSI 

J 2 u 

SSI, SS2 b , SS2# 

J-i 

SSI, SS2 b , SS2* 1 , SS3 


Table 4 Existence of steady-states in the regions of the op¬ 
erating diagrams of Fig. [3|i) and Fig. [ 6 ])i). 




Fig. 4 The curves A (black), A (red) and A (green) for 
case (b). (i) : regions of steady-state existence, with mainte¬ 
nance. (ii) : regions of steady-state existence and their stabil¬ 
ity, without maintenance. On the right, a magnification for 
0 < D < 0.1 


Region 

Steady states 

Ji 

SSI 

Ji 

SSI, SS2 b , SS2# 

J 3 

SSI, SS2 b , SS2#, SS3 


Table 5 Existence of steady-states in the regions of the op¬ 
erating diagram of Fig[4])i). 


4.2 Operating diagram: case (b) 

This case corresponds to the parameter values used 
by nm, except that Ks,h 0 , c is changed from 1.0 x 10 6 
to 4.0 x 10 -6 . We have seen in Table[2]that the curves 
.Ti and r 2 are defined for A < A 1 and A < A 2 , re¬ 
spectively and A (A) < F 2 {D) for all D < A 2 , where 
D\ = 0.329 and D 2 = 0.236. Therefore, they sepa¬ 
rate the operating plane (^ch.in, D) in three regions, as 
shown in Fig. |4^i) , labelled J \, J 3 and 

The results are summarised in Table [5j which shows 
the existence of the steady-states SSI, SS2 and SS3 in 
the regions of the operating diagram in Fig. |4])i). Note 
that the region J 2 has disappeared. 
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Fig. 5 The curve Ti for case (c). (i) : regions of steady-state 
existence, with maintenance, (ii) : regions of steady-state ex¬ 
istence. without maintenance and their stability. On the right, 
a magnification for 0 < D < 0.1. 


Region 

Steady states 

Ji 

SSI 

Ja 

SSI, SS2 b , SS2* 1 


Table 6 Existence of steady-states in the regions of the op¬ 
erating diagram of Fig [5ji). 


4.3 Operating diagram: case (c) 

This case corresponds to the parameter values used 
by m, except that K,s,h 2 .c is changed from 1.0 x 10 6 
to 7.0 x 10 -6 . We have seen in Table [ 2 ] that the curve 
A is defined for D < D\ = 0.287 and that I 2 is empty 
so that SS3 does not exist. Therefore, A separates the 
operating plane (SAinA) in two regions, as shown in 
Fig. [5ji) , labelled A and A- 

The results are summarised in Table [6j which shows 
the existence of the steady-states SSI and SS2 in the 
regions of the operating diagram in Fig. [5ji). Note that 
the region A of existence of SS3 has disappeared. 




Fig. 6 The curves i~i (black), I2 (red) and T3 (green) for 
case (d). (i) : regions of steady-state existence, with mainte¬ 
nance. (ii) : regions of steady-state existence and their stabil¬ 
ity, without maintenance. On the right, a magnification for 
0 < D < 0.1. 

that they are tangent for D = Z) 3 , where D 3 = 0.238, 
-^2 min = 0.101, D 2m ax = 0.198 and D 3 = 0.161. There¬ 
fore, A and A separate the operating plane (Ah,in, D) 
in four regions, as shown in Fig. [6](i) , labelled A , A ■ A 
and A- The results are summarised in Table [4j which 
shows the existence of the steady-states SSI, SS2 and 
SS3 in the regions of the operating diagram in Fig.[6](i). 


4.5 Stability of steady-states 

We know that SSI is always stable. The analytical study 
of the stability of SS2 and SS3 is very difficult because 
the conditions for Routh-Hurwitz in the 6-dimensional 
case are intractable. For this reason we will consider in 
Section [5] the question of the stability only in the case 
without maintenance, since the system reduces to a 3- 
dimensional. The general case will be considered only 
numerically in Section [6] 


4.4 Operating diagram: case (d) 

We end this discussion on the role of kinetic param¬ 
eters by the presentation of this case, which presents 
a new behaviour that did not occur in the preceding 
cases: there exists D 2 min such that for D < D 2 min the 
system cannot have a positive steady-state SS3. This 
case corresponds to the parameter values used by m, 
except that three of them are changed as indicated in 
Table [3] This table shows that the curves A and A 
are defined for D < D\ and D 2m in < D < D 2max and 


5 Local stability without maintenance 

When maintenance is not considered in the model, the 
steady-states SSI, SS2 and SS3 are given by 

1. SSI = (0,0,0,4 n ,0,0) 

2. SS2 = Ao, xi, 0, so, Si, s 2 ) where S 2 a solution of 
equation 


So" =i>{s 2 ) 


Mq(D , s 2 ) + 


Mi(D, s 2 ) + s 2 

1 ~ id 
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and 


So = Mq(D, s 2 ), Si = Mi(D, S 2 ) 

xq = s™ - s 0 , X\ = s' 0 n - s 0 - si (47) 

3. SS3 = (xo, Xi, X 2 , so, si, S 2 ) where 

S 2 = M 2 (D), sq = M 0 (D,s 2 ), si = Afi(D, S 2 ) 
x 0 = s™ - s 0 , x 1 = so 11 - s 0 - si (48) 

X2 = (1 - w) (s™ - S 0 ) - Si - S 2 

Proposition 1 Let SS2 = (xo, Xi, 0, so, Si, S 2 ) be a 
steady-state. Then SS2 is stable if, and only if, p 2 (s 2 ) < 
D and 4^- > 0. 

as 2 

Therefore, SS2 b is always unstable and SS2 # is stable if, 
and only if, ia 2 (s 2 ) < D. This last condition is equiv¬ 
alent to M 2 (D) > s\, which implies that F 3 (D) > 0. 
Hence, if SS3 exists then SS2 # is necessarily unstable. 
Therefore, SS2® is stable if, and only if, F 3 (D) > 0 and 
SS3 does not exist. 

Proposition 2 Let SS3 = (xo, xi, X 2 , so, si, s 2 ) be a 
steady-state. If F 3 (D) > 0 then SS3 is stable as long 
as it exists. If F 3 (D) < 0 then SS3 can be unstable. 
The instability of SS3 occurs in particular when s 2 *s 
sufficiently close to s b 2 , that is to say SS3 is sufficiently 
close to SS2 11 . 


The condition F 3 (D) > 0 is equivalent to 

jfc{M 2 (D)) > 0, that is to say s 2 = M 2 (D) £ [s^Sj). 

If Ts 2 < 0, that is to say s 2 £ (s 2 ,s 2 ), then SS3 can be 
unstable. 

When D is such that F 3 (D) < 0, the determination 
of the boundary between the regions of stability and 
instability of SS3 needs to examine the Routh-Hurwitz 
condition of stability for SS3. For this purpose we define 
the following functions. Let SS3 = (xo, Xi, x 2 , sq, s 1, s 2 ) 
be a steady-state. Let 


p, _ Qpp Q _ 

9sn ’ dsn ’ 0si ! 


_ dPl 


H=- 


1 = 


_ d F -2 


evaluated at the steady-state SS3 defined by (481, that 
is to say, for 


s 2 — M 2 (D), so — Mq(D , s 2 ), s 1 -M 1 {D,s 2 ) 

For D £ I 3 and s™ > F 2 {D), we define 

F 4 (D, 4 n ) = {EIxox 2 + [E{G + H) - (1 - lj)FG] x 0 xi)/ 2 

+ (Ix 2 + (G + H)xi + ujFxq)GIxiX 2 (49) 


where f 2 = Ix 2 -\-(G+H)xi + (E+(jjF)xq. Notice that to 
compute F 4 ( D , s™), we must replace xo, Xi, x 2 , So, si 
and s 2 by their values at SS3, given by (48). Hence, this 
function depends on the operating parameters D and 



Existence 

Stability 

SSI 

Always exists 

Always stable 

SS2 b 

>F 1 (D) 

Always unstable 

SS2 # 

s i 0 n >F 1 (D ) 

F 3 (D) > 0 and s" 1 < F 2 (D) 

SS3 


F 3 (D) > 0 or 

F 3 (D) < 0 and F 4 (D, 4 n ) > 0 


Table 7 Existence (with or without maintenance) and sta¬ 
bility (without maintenance) of steady-states. 


s), n . For each fixed D £ I 3 , F (£>, Sq 11 ) is polynomial in 
s), n of degree 3 and tends to +00 when s™ tends to +00. 
Therefore, it is necessarily positive for large enough s™. 
The values of the operating parameters D and Sq 11 for 
which F (D, Sg n ) is positive correspond to the stability 
of SS3 as shown in the following proposition. 

Propositions Let SS3 = (xo, Xi, x 2 , So, Si, s 2 ) be a 
steady-state. If F 3 (D) < 0 then SS3 is stable if, and 
only if, F 4 (As™) > °- 

The results on the existence of steady states (with 
or without maintenance) of Lemma |TJ Lemma [ 2 ] and 
Lemma [ 3 } and their stability (without maintenance) of 
Prop|T] Prop [2] and Prop[3j are summarised in Table [7] 


5.1 Operating diagram: case (a) 

This case corresponds to the parameter values used 
by [10. but without maintenance. We see from Table 
[ 2 ] that the curves r 3 and F 2 of the operating diagram, 
given by Eq. [45] and Eq. [46j respectively, are defined 
now for D < Di = 0.452 and D < D 2 = 0.393, respec¬ 
tively and that they are tangent for D = D 3 = 0.078. 
Beside these curves, we plot also on the operating dia¬ 
gram of Fig.[3jii), the curve J 3 of equation 

F 4 (D,Y 3 Y4S ch - m )=0 (50) 

According to Prop. [3] this curve is defined for D < 
D 3 = 0.078 and it separates the region of existence 
of SS3 into two subregions labelled J 3 and J 3 , such 
that SS3 is stable in J 3 and unstable in J 3 . The other 
regions J\, J 2 and J 4 are defined as in the previous 
section. The operating diagram is shown Fig. |gii). It 
looks very similar to Fig.[3])i), except near the origin, as 
it is indicated in the magnification for 0 < D < D 3 = 
0.078. From Table [7J we deduce the following result 

Proposition 4 Table^ shows the existence and stabil¬ 
ity of the steady-states SSI, SS2 and SS3 in the regions 
of the operating diagram in Fig. [3|(n,). 
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Region 

SSI 

SS2 b 

SS2# 

SS3 

Ji 

s 




Ji 

s 

U 

S 


J:i 

s 

U 

u 

s 

Ja 

s 

u 

u 


J 5 

s 

u 

u 

u 


Table 8 Existence and stability of steady-states in the re¬ 
gions of the operating diagrams of Fig. iii) and Fig. [efii). 


Region 

SSI 

SS2 b 

SS2» 

SS3 

Ji 

S 




J 3 

S 

U 

U 

s 

Ja 

S 

U 

u 


Js 

S 

U 

u 

u 


Table 9 Existence and stability of steady-states in the re¬ 
gions of the operating diagram of Fig. 


Region 

SSI 

SS2 b 

SS2» SS3 

Ji 

S 



Ja 

S 

U 

U 


Table 10 Existence and stability of steady-states in the re¬ 
gions of the operating diagram of Fig. 0iO- 


5.2 Operating diagram: case (b) 


We see from Table [2] that the curves A and A are de¬ 
fined now for D < Di = 0.349 and D < Z) 2 = 0.256, 
respectively and that F\(D) < A(-D) for all D. Beside 
these curves, we plot also on the operating diagram of 
Fig.0ii), the curve r 3 of equation (Eq. 501 which sepa¬ 
rates the region of existence of SS3 into two subregions 
labelled J 3 and J 3 , such that SS3 is stable in J 3 and 
unstable in J 3 . Therefore, the curves A, A and A sep¬ 
arate the operating plane (S' c h,m,-D) into four regions, 
as shown in Fig. 0ii), labelled J \, J 3 , J,\ and J 5 . 


Proposition 5 Table^shows the existence and stabil¬ 
ity of the steady-states SSI, SS2 and SS3 in the regions 
of the operating diagram in Fig. ^ii) 


5.3 Operating diagram: case (c) 

We see from Table [2] that A is defined for D < D\ = 
0.303 and that I 2 is empty so that SS3 does not exist. 
Therefore, A separates the operating plane (Ah,in,-A 
into two regions, as shown in Fig. 0ii), labelled J\ and 
Ja- 

Proposition 6 Table \10\ shows the existence and sta¬ 
bility of the steady-states SSI, SS2 and SS3 in the re¬ 
gions of the operating diagram in Fig. 


5.4 Operating diagram: case (d) 


We see in Table [3] that the curves A and A are de¬ 
fined for D < Di and A 2m j„ < D < Z3 2max and that 
they are tangent for D = D 3 , where D i = 0.258 and 
F>imin = 0.121, D^max = 0.218 and D 3 = 0.181. Beside 
these curves, we plot also on the operating diagram of 
Fig. 0ii), the curve A defined by Eq. 50 which sepa¬ 
rates the region of existence of SS3 into two subregions 
labelled J 3 and J 3 , such that SS3 is stable in J 3 and 
unstable in J 3 . Therefore, the curves A, A and A sep¬ 
arate the operating plane (S c h,i n ,-D) into five regions, 
as shown in Fig. 0ii), labelled J\, Ji, J 3 , Ja and J 3 . 


Proposition 7 Table^ shows the existence and stabil¬ 
ity of the steady-states SSI, SS2 and SS3 in the regions 
of the operating diagram in Fig. [6^ ii). 


6 Numerical analysis to confirm and extend the 
analytical results 

The aim of this section is to study numerically (the 
method is explained in Appendix [0) the existence and 
stability of the steady-states SS2 and SS3. We obtain 
numerically the operating diagrams that were described 
in Sections 0 and 0 The results in this section confirm 
the results on existence of the steady-states obtained 
in Section 0 in the case with or without maintenance 
and the results of stability obtained in Section [5] in the 
case without maintenance. These results permit also to 
elucidate the problem of the focal stability of SS2 and 
SS3, which was left open in Section [475} 


6.1 Operating diagram: case (a) 

We endeavoured to find numerically the operating con¬ 
ditions under which SS3 is unstable, previously unre¬ 
ported by pTO] ■ Given that we have determined analyti¬ 
cally in Proposition [ 2 ] that when SS3 is close to SS2^ it 
becomes unstable, we performed numerical simulations 
with the parameters defined in Table 0 over an oper¬ 
ating region similar to that shown in Fig. 2 from m 
whilst also satisfying our conditions. In Fig. 0 we show 
the case when maintenance is excluded. When magni¬ 
fied, we observe more clearly that region J 5 does exist 
for the conditions described above, and also note that 
the region Ja occurs in a small area between J\ and J 3 , 
which corresponds to the results shown in Fig. 3[ii), 
and is in agreement with Proposition 0 In Fig. 8j we 
confirm that region J 3 does exist for the conditions 
described above, when maintenance is included, but 
could not be determined analytically, the curve A is 
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absent in Fig. |3ji) . Furthermore, we demonstrate that 
a Hopf bifurcation occurs along the boundary of F 3 (D) 
for values of D < D 3 by selecting values of 5 c h,i n (indi¬ 
cated by (a) — (<5) in Fig. [8]) at a fixed dilution rate of 
0.01 d _1 , and running dynamic simulations for 10000 d. 
The three-dimensional phase plots, with the axes rep¬ 
resenting biomass concentrations, are shown in Fig. [9j 
and show that as S C h,in approaches J 3 from J 3 , emer¬ 
gent periodic orbits are shown to diminish to a stable 
limit cycle at the boundary (see Appendix [Bjfor proof). 
Subsequently, increasing S' c h,in to J 3 results in the orbit 
reducing to a fixed point equilibrium at SS3. 



Fig. 7 Numerical analysis for the existence and stability of 
steady-states for case (a), without maintenance. On the right, 
a magnification for 0 < D < 0.16. 


(«) ( 0 ) 




Fig. 9 Three-dimensional phase plane diagrams of the 
biomass dynamics for t = 10000 d, showing initial (green dot) 
and final (red dot) conditions for a dilution rate, D = 0.01 d -1 
and chlorophenol input, S c h,in ( kgCOD/m 3 ) of a) 0.01 - the 
system converges to SSI, /3) 0.097 - the system enters a peri¬ 
odic orbit of increasing amplitude, ultimately converging to 
SSI, 7 ) 0.10052 - the system is close to a stable limit cy¬ 
cle, 5) 0.16 - the system undergoes damped oscillations and 
converges to SS3. 


6.2 Operating diagram: case (b) 



Fig. 8 Numerical analysis for the existence and stability of 
steady-states for case (a), with maintenance. This is a mag¬ 
nification for 0 < D < 0.1, showing the presence and extent 
of region 77s undetectable by the analytical method. The co¬ 
ordinates labelled (a) — (5) are subsequently used to simulate 
the system dynamics, as shown in the proceeding Fig. [9] 


Whilst the numerical parameters chosen for this work 
are taken from the original study [TO] , there some¬ 
what arbitrary nature leaves room to explore the im¬ 
pact of the parameters on the existence and stability 


of the steady-states. Case (b), discussed in Sections 4.2 


and |5.2[ involves a small increase to the half-saturation 
constant (or inverse of substrate affinity), -Ksh 2iC , of 
the chlorophenol degrader on hydrogen. Following the 
same approach as with the preceding case, we confirm 
in Fig. 10 4) the Proposition [ 5 ] in the scenario without 
maintenance. Furthermore, the extension of this propo¬ 
sition with maintenance included, corresponding to the 
existence and stability of all three steady-states given 
in Table |9j is show in Fig |To](ii) . It shows the region 
J 5 that cannot be obtained analytically (cf. Fig. Si))- 
In both cases, region J 2 has disappeared, as observed 
analytically. Additionally, the ideal condition J 3 , where 
all organisms are present and stable, diminishes. 


6.3 Operating diagram: case (c) 

Here, Ks fj 2 ,ci was further increased and confirm the 
Proposition [6j where the function SS3 never exist and 
SS2 never stable for the case without maintenance. The 
extension of this proposition to the case with mainte- 
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nance, shown in Table 10 produce similar results as 
shown in the comparison of Figs. 11 i) and (ii). 


6.4 Operating diagram: case (d) 

With the final investigated scenario, where k m ^ 2 < 
km ,ch and A(s,h 2 < -^s,h 2 ,c> we observe once again the 
presence of all operating regions, J\ — without and 
with maintenance, as shown in Fig. |12| It can be seen 
that regions Ja and increase at low dilution rates 
across a much larger range of S'ch.in than in the default 
case (a), and the desired condition (stable SS3) is re¬ 
stricted to a much narrower set of D. 

As with the previous cases, the numerical analysis 
for case (d) confirms the Proposition [7] without mainte¬ 
nance and its extension to the case with maintenance, 
indicated in Table |U 



Fig. 10 Numerical analysis for the existence and stability of 
steady-states for case (b). (i) : without maintenance, (ii) : 
with maintenance. 


7 The role of kinetic parameters 

Finally, we give brief consideration to the characteri¬ 
sation of the four cases discussed in the preceding sec¬ 
tions. The main difference between cases (a) or (b) and 
cases (c) or (d) is that, for small values of D, the co¬ 
existence steady-state SS3 can exist for cases (a) and 
(b), but cannot exist for cases (c) or (d). The cases 
(a) or (b) occur if and only if s 2 (0) < M 2 ( 0) holds or 
s°(0) = M 2 (0) and (0) < ^=y(0) hold, that is to say 

Lqclo R 2 a 2 

- < - or 

TO 0 — a 0 m 2 — a 2 

(51) 

L 0 ao _ K 2 a 2 L 0 m 0 K 2 m 2 

Too - a 0 m 2 - a 2 (m 0 - a 0 ) 2 ( m 2 - a 2 ) 2 

(52) 



Fig. 11 Numerical analysis for the existence and stability of 
steady-states for case (c). (i) : without maintenance, (ii) : with 
maintenance. On the right, a magnification for 0 < D < 0.1. 



S'cli, i ii <Sch,in 

Fig. 12 Numerical analysis for the existence and stability of 
steady-states for case (d). (i) : without maintenance, (ii) : 
with maintenance. 


The cases (c) or (d) occur if and only if s 2 (0) > M 2 (0) 
holds or s°(0) = M 2 (0) and ^f(0) > ^yf(O) hold, that 
is to say 


Lpap 
to 0 - a 0 


Lqciq A 2 a 2 

> - or 


to 0 - a 0 m 2 — a 2 


(53) 


K 2 m 2 

m 2 - a 2 (?7i 0 - a 0 ) 2 " (m 2 - a 2 ) 


K 2 a 2 L 0 m 0 

and 


> 


(54) 


Notice that it is easy to make the difference between 
case (c) and case (d): the first occurs when M 2 (Di) < 
s 2 (D i) and the second when M 2 (D{) > s^-Di). Since 
D\ is the positive solution of the algebraic quadratic 
equation s 2 (H) = s 2 (Zl), it is possible to have an ex¬ 
pression for D\ with respect to the biological param- 
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eters. However, this is a complicated expression in¬ 
volving many parameters and the preceding conditions 
M 2 (.Di) < s 2 (-Oi) or > s 2 (Di) have no bio¬ 

logical interpretation. We simply remark here that the 
function s^(D) has a vertical asymptote for D = m 0 —aQ 
and the function M 2 (D) has a vertical asymptote for 
D = m 2 — a 2 . Therefore, if Too — ao < m 2 — a 2 then 
case (c) occurs, so that a necessary (but not sufficient) 
condition for case (d) to occur is mo — a 3 > m 2 — a 2 ■ If 
m 2 is sufficiently small then case (d) can occur. 

The observations from the numerical analysis sug¬ 
gest that the role of the chlorophenol degrader as a 
secondary hydrogen scavenger is critical in maintain¬ 
ing full chlorophenol mineralisation and system stabil¬ 
ity, particularly at higher dilution rates, as shown by 
comparing cases (c) and (d) . More significantly, the 
results coupled with the parameter relationships shown 
in Eqs. 5ffl54 highlight the necessary conditions under 
which the ideal case (SS 3 stable) is achieved and, in 
general, this is a coupling of the two key parameters 
describing the half-saturation constant and maximum 
specific growth rates between the two hydrogen com¬ 
petitors. 


8 Conclusions 

In this work we have generalised a simplified mecha¬ 
nistic model describing the anaerobic mineralisation of 
chlorophenol in a two-step food-web. We are able to 
show complete analytical solutions describing the exis¬ 
tence and stability of the steady-states in the case that 
maintenance is excluded from the system, whilst with 
a decay term present, purely analytical determination 
of stability is not possible. 

We confirm the findings of previous numerical anal¬ 
ysis by m that with chlorophenol as the sole input 
substrate, three steady-states are possible. However, 
the analysis goes further and we determine that under 
certain operating conditions, two of these steady-states 
(SS2 and SS3) can become stable, whilst SSI always 
exists and is always stable. Furthermore, without main¬ 
tenance we can explicitly determine the stability of the 
system, and form analytical expressions of the bound¬ 
aries between the different stability regions. 

As the boundary of J 3 is not open to analytical 
determination in the case with maintenance, we de¬ 
termined numerically (substituting the general growth 
function with the classical Monod-type growth kinet¬ 
ics) the existence and stability of the system over a 
range of practical operating conditions (dilution rate 
and chlorophenol input). For comparison and confirma¬ 
tion, we also performed this for the case without main¬ 
tenance and found the same regions in both cases, with 


variations only in their shape and extent. For example, 
whilst the boundary between J\ and 3 4 terminates at 
the origin without maintenance, with maintenance it 
is located at ^(O)/!^^ « 0.0195. More interestingly, 
the addition of a decay term results in an extension 
of the SS3 unstable steady-state, reducing the poten¬ 
tial for successful chlorophenol demineralisation at rel¬ 
atively low dilution rates and substrate input concen¬ 
trations. Additionally, we show that at the boundary 
between J 3 and J7s, a Hopf bifurcation occurs and a 
limit cycle in SS3 emerges. 

Finally, we gave an example of how the model could 
be used to probe the system to answer specific ques¬ 
tions regarding model parameterisation. Here we have 
indicated that a switch in dominance between two or¬ 
ganisms competing for hydrogen results in the system 
becoming unstable and a loss in viability. This is per¬ 
haps intuitive to microbiologists, but here it has been 
proven using mathematical analysis, and could be used 
to determine critical limits of the theoretical parameter 
values in shifting between a stable and unstable system. 
Whilst parameters are not arbitrary in real organisms, 
the potential for microbial engineering or synthetic bi¬ 
ology to manipulate the properties of organisms makes 
this observation all the more pertinent. 


A Numerical methods 

We consider sets of operating parameters (D and 5 c h,i n ) f° r 
each of the three steady-states, and using Matlab, the com¬ 
plex polynomials for each steady-state can be solved by sub¬ 
stitution of parameter values (see Table [lj into the explicit 
solution. By investigating the signs of the solutions and the 
eigenvalues, respectively, we determine which steady-states 
are meaningful and stable. By exploring a localised region 
of suitable operating parameters, we then generate a phase 
plot showing where each steady-state is stable, bistable or 
unstable. 


B Proof for Hopf Bifurcation 

In Section [ 6 ] we show the operating diagrams with the pa¬ 
rameters given in Table [T] and determine numerically that 
as the parameter S c h,i n increases at a fixed dilution rate 
(D = 0.01 d _1 ), the system bifurcates through several sta¬ 
bility domains. We claim that as we cross the boundary be¬ 
tween regions J 5 and J 3 , we observe a Hopf bifurcation, and, 
in J$, close to the boundary with a limit cycle appears. 
In order to test this numerically, we checked the real parts of 
the six eigenvalues at each point along the transect shown in 
Fig[8](10000 points in total), and plotted their values. Fig |13| 
indicates the conditions for a Hopf bifurcation are satisfied 
as eigenvalues 2 and 3 both change their sign when passing 
through the coordinate (0, 0.1034) and the real part of all 
eigenvalues 1, 4 and 6 remain negative. 




Eigenvalue 5 Eigenvalue 3 Eigenvalue 1 
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Fig. 13 Real parts of the eigenvalues determined at D = 0.01 
and 5 c h,in = [0.08,0.12], in the case with maintenance. The 
red vertical lines indicate the location where the eigenvalue 
crosses zero. 


C General case 


As mentioned at the end of Section [2] our study does not re¬ 
quire that growth functions are of Monod type (Eq. |23| ). Ac¬ 
tually, the results are valid for a more general class of growth 
functions satisfying the following conditions, which concur 
with those given by Eq. |23| 

HI For all so > 0 and S 2 > 0 then 0 < mo (so, s 2 ) < +00 and 
/.to (0, s 2 ) = 0, /ro (so, 0) = 0. 

H2 For all si > 0 and S 2 > 0 then 0 < fii (s 1 , S 2 ) < +00 and 

Mi (0, S 2 ) = 0. 

H3 For all S 2 > 0 then 0 < M 2 (• 52 ) < +00 and y. 2 (0) = 0. 

H4 For all so > 0 and S 2 > 0, 


d/10 

dso 


(so, S 2 ) > 0, 


<3mo 

ds 2 


(so,S 2 ) > 0. 


H5 For all si > 0 and S 2 > 0, 


9/j.i 
ds 1 


(si,s 2 ) > 0, 


<9mi 

ds 2 


{si,s 2 ) < 0. 


H6 For all s 2 > 0, (S 2 ) > 0. 

ds 2 

H7 The function S 2 1 —>■ /1 0 (+ 00 , 52 ) is monotonically increas¬ 
ing and the function S 2 m- mi(+°°> s 2 ) is monotonically 
decreasing. 

We use Eq. |31| Eq. |32| and Eq. |33| to define Mo(y,s 2 ), 
Mi(y,s 2 ) and M 2 (y), respectively. 


Lemma 4 Let s 2 > 0 be fixed. There exists a unique function 


y e [0, mo(+°°, s 2 )) M 0 (y , s 2 ) e [0, + 00 ), 

such that for so > 0, S 2 > 0 and y £ [0, mo(+oo, s 2 )), we have 
so = M 0 (y, s 2 ) </=> y = yo(s 0 ,s 2 ) (55) 


Lemma 5 Let s 2 > 0 be fixed. There exists a unique function 

y £ [ 0 ,mi(+°o,S 2 )) s-t Mi (y, s 2 ) £ [ 0 , +00), 

such that for si > 0, S2 > 0 and y £ [0, £ [0, pi (+ 00 , S 2 )), we 
have 


si = Mi(y,s 2 ) <!=> y = yi(si,s 2 ) 


(56) 


Lemma 6 There exists a unique function 

y 6 [0, M 2 (+oc)) s-t M 2 (y) £ [0, +oo), 
such that, for s 2 > 0 and y £ [0, M2(+oo)) we have 
S2 = M 2 (y) •<=> y = M2(s2) 


(57) 


We use Eq. 35 to define the functions Sj(D) and s\(D) 

Lemma 7 For D + ao < Mo(+ 00 ,+ 00 ) and D a\ < 
Ml(+oo,0) there exist unique values stj and s\ such that 


Mo(+°o, s°) = D + ao, y,i(+oo,sl) = D + ai 


(58) 


Let ui < 1. We use Eq. |34| to define i/i(s 2 , D) in the general 
case: we let : (s§, Sj) —» R defined by 


ip(s 2 ) = Mo(D + ao, s 2 ) + 


Mi (D + ai,s 2 ) + S2 

1 — to 


(59) 


It should be noted that i/;(s 2 ) > 0 for s§ < S 2 < Sj. From 
Eq. |55| Eq. |56| and Eq. |58| we deduce that 

Mo(D + ao, s®) = + 00 , Mi(D + ai, S2) =+00 

Therefore, we have 

lim ip(s 2 ) = lim if(s 2 ) = +00 
s 2 ->s§ s 2 ->s 2 


Hence, the function i/)(s 2 ), which is positive and tends to +00 
at the extremities of the interval (s^Sj), has a minimum 
value on this interval. We add the following assumption: 

H8 The function if has a unique minimum s 2 on the interval 
(sf^Sj) and ^^{s 2 ) is negative on (sj,S 2 ) and positive on 
(s 2 , S 2 ), respectively. 

The function ip together with the values Sj, Sj and s 2 
all depend on D. However, to avoid cumbersome notations 
we will use the more precise notations ip(s 2 ,D), s^D), s\(D) 
and S 2 (D) only if necessary. 

We use Eq. |36| Eq. |38| and Eq. |39| to define Fi(D), F 2 (D) 
and Fs(D) in the general case: 


Fi (D) = inf ip(s 2 ) = ip (s 2 ) 

(60) 

s 2 e(s°,si) 


F 2 (D) = iP(M 2 (D + a 2 )) 

(61) 

F 3 (D) = ^ (M 2 {D + a 2 )) 
as 2 

(62) 


The function Fi ( D ) is defined for 

D £ h = {D > 0 : s%(D) < 4(D)} 


The function F 2 (D) and Fo(D) are defined for 

D £ I 2 = {D £ h : sg(D) < M 2 (D + a 2 ) < sl{D)} 

For all for D £ I 2 , Fi(D) < F 2 (D). The equality Fi(D) = 
F 2 {D ) holds if, and only if, M 2 (D + a 2 ) = s 2 (D) that is, 
( M 2 (D + a 2 )) = 0, that is if, and only if, Fs(D) = 0. 

As it will be shown in Appendix [DJ the Lemmas [ 2 ] and 
[3j stated in Section [5] in the particular case of the Monod 
type growth functions (Eq. |23| ) , are true in the general case 
of growth functions satisfying assumptions H1-H8. 
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D Proofs 

In this Section we give the proofs of the results. In these 
proofs, we do not assume that the growth function are of 
Monod type (Eq. |23| ). We only assume that the growth func¬ 
tions satisfy H1-H8. 


D.l Existence of steady-states 


Proof [Lemma [3] Since xo > 0, x\ > 0 and X 2 > 0, then, as a 
consequence of Eq. |25| Eq. |26| ) and Eq. |27| we have 

Mo(so, S 2 ) = D + a, Ml( s l, s 2 ) = D + b, ^ 2 (^ 2 ) = D + c 

Hence, so, si and S 2 are given by Eq. |43| Using Eq. |28| Eq. |29| 
and Eq. |30| we have Eq. |44| For X 2 to be positive it is necessary 
that so, si and S 2 satisfy the condition 

(1 - to)(sQ n - so) > si + S 2 , (65) 


Proof [Lemm a [T] Assume first that xq = 0. Then, as a conse¬ 
quence of Eq. |28| we have so = sJf 1 and, as a consequence of 
Eq. |29[ we have 

Ds 1 + fii(si, S 2 )xi = 0 

which implies si = 0 and Ml ( s i, s 2 )tri = 0. Therefore, as a 
consequence of Eq. |26| we have xi = 0. Replacing xo = 0 and 
x\ = 0 in Eq. |30| we have 

DS2 + fl2(s2)X2 = 0 


which implies S 2 = 0 and M 2 (s 2 )x 2 = 0. Therefore, as a con¬ 
sequence of Eq. |27| we have X 2 = 0. Hence, the steady-state 
is SSI. 

Assume now that xi = 0. Then, as a consequence of Eq. |30[ 
we have 

Ds 2 + aj/i 0 (so,S 2 )rro + M 2 (s 2 )x 2 = 0 

which implies S2 = 0 , mo(so, S 2)xo = 0 and M2(s2)£2 = 0 . 
Therefore, as a consequence of Eq. |25| we have xq = 0. As 
shown previously this implies that the steady-state is SSI. 
Evaluated at SSI the Jacobian matrix of Eqs. Em is 


‘ — D — ao 0 
0 — D — ai 

0 0 

0 0 

0 0 

0 0 


0 000- 

0 0 0 0 

—D — 02 0 0 0 

0 -D 0 0 

0 0 -D 0 

0 0 0 —D_ 


Thus, SSI is stable. □ 

Proof [Lemm a [2] Since xo > 0 and x\ > 0, then, as a conse¬ 
quence of Eq. |25| and Eq. |26| we have 


Mo(so,S 2 ) = D + a 0 , A t 1 (si,s 2 ) = D + ai 


Hence, we have 


so = Mo(D + ao, S 2 ), si = Mi(D + ai, S 2 ) (63) 


Using Eq. |28| and Eq. |29| we have Eq. |42| Using Eq. |30| we 
have 


—32 + (4 n - SO - Si) - Ul(s o' - so) = 0 


(64) 


If to > 1 this equation has no solution. If to < 1 this equation 
is equivalent to 

in _ . SI + S2 

s 0 — so H— -■ 

1 — U! 

Using Eq. |63| we see that S 2 must be a solution of Eq. |41| Since 
si > 0 and S 2 > 0 then, from Eq. |64| we have necessarily 

si + s 2 = (1 - d))(so n - so) > 0 


so that Sp n — so > 0. From Eq. |42| 
s q 1 — sq > 0 and S 2 > 0 then, from Eq. 


we de duc e that xo > 0. Since 
we have necessarily 


64 


w(so 1 - s 0 ) + s 2 = sg 1 - s 0 - si > 0 


so that si 11 — sq — si >0 From Eq. 


42 


we deduce that xi > 0 . 

□ 


If cj > 1 this equation has no solution. If ui < 1 this equation 
is equivalent to the condition 


SfJ 1 > so T 


Sl + S2 


1 — to 

Using Eq. |43| this condition is the same as 

sj, n >4>(M 2 (D + a 2 )) = F 2 (D) 


Therefore, from Eq. 65 we have Sq 11 —so > 0 and Sg 1 — so — si > 


0 , so that xq > 0 andxi > 0 . 


D.2 Stability of steady-states 


We use the change of variables 


20 = SO +3:0, 2l = Sl+Xl—X 0 , Z 2 = S2+X2+0JX0-X! 

Therefore, Eqs. Em with ao = ai = a ,2 = 0 , become 
dx 0 


= —Dxg + fl 0 (zo - X 0 ,Z 2 - OJX 0 + Xl - X 2 ) X 0 
= —Dx 1 + /ll (2:1 + x 0 — Xi, 22 - tox 0 + Xl — x 2 ) Xi 
= —Dx 2 + M 2 (z 2 — COXo + Xl - X 2 ) X2 
= D - zo) 


d t 

dxi 

d t 
dx 2 
d t 
dz 0 
dt 
dzi 

-aT = 


dz 2 

dt 


= —DZ2 


( 66 ) 

(67) 

( 68 ) 

(69) 

(70) 

(71) 

(72) 


In the variables (xo, xi, X2,20, zi, 22) where 21, 2:2 and 23 are 
defined by Eq. | 66 | the steady-states SSI, SS2 and SS3 are 
given by 

1. SSI = (0, 0, 0, sb n , 0,0) 

2. SS2 = (xo, xi, 0, Sq 11 , 0, 0), where xo and xi are defined by 

Eq-El 

3. SS3 = (xo, xi, X2, Sq 11 , 0, 0), where xo, xi and X2 are de¬ 
fined by Eq. |48| 

Let (xq, xi , X2, sy , 0 , 0 ) be a steady-state. The Jacobian 
matrix of Eqs. |67|72| has the block triangular form 

t _ [Jr J 2 ' 

0 J 3 


where 


Ji 


Mo - D - (E + ljF)x 0 Fx 0 

(G + uiH)x! mi - D - (G + ff)n 
— Ldlx 2 lx 2 


— Fx 0 
Hx 1 

H 2 ~ D — 1 X 2 


Ex 0 

0 

Fx 0 


—D 

0 

0 

0 

Gxi 

-Hx 1 

, J 3 = 

0 

—D 

0 

0 

0 

Ix 2 


0 

0 

—D 
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and 

E= d }f £ f= 9 ^- G=— H =-—- /= — 

dso ’ ds 2 ' ds 4 ’ &S 2 ’ ds 2 

are evaluated at the steady-state. 

Since J is a block triangular matrix, its eigenvalues are 
— D (with multiplicity 3) together with the eigenvalues of the 
3x3 upper-left matrix Ji. Note that we have used the oppo¬ 
site sign for the partial derivative H = —dfi 4 /ds 2 , so that all 
constants involved in the computation become positive, which 
will simplify the analysis of the characteristic polynomial of 
Jr. 


Proof [Proposition [I] Evaluated at SS2, the matrix Jr is 


Jr 


— (Ef-ujF)x 0 Fx 0 —Fx 0 

(G + ujH) xt -(G + H)x 1 Hx 1 
0 0 fi 2 — D 


Since Ji is a block triangular matrix, its eigenvalues are sim¬ 
ply ii 2 —D, together with the eigenvalues of the 2x2 upper-left 
matrix. Note that the trace of this 2x2 matrix is negative. 
Hence, its eigenvalues are of negative real part if, and only if, 
its determinant is positive, that is if, and only if, 


E(G + H) - (1 - lo)FG > 0 (73) 

Using 

dMp 
ds 2 

dMi 
ds 2 

we deduce from 

U \ a, ^ , M\(D,S2 ) + S2 

’4>{s2) = Mq(D, s 2 ) H--- 

1 — ui 

that 

dip _ dMp | + 1 _ F ! g + 1 

dS2 dS2 1 — LJ E 1 — (Jj 

Hence, 


dpo 

ds2 

_ dpi 

ds2 


dfio 


dso 

dpi 

dsp 


= —F/E 


= H/G 


d± = E(G + H) - (1 - uj)FG 

ds 2 (1 -u)EG 1 ' 

Therefore, the condition of stability, (Eq. |73| ), is equivalent 
to > 0. Hence, we have proved that SS2 is stable if, and 

only if, ^ 2 ( 32 ) < D and ^ >0. □ 


Proof [Proposition [ 2 ] Evaluated at SS3, the matrix Ji is 


Ji 


— (E + uF)xo Fx 0 —Fx 0 

(G + ujH)x r ~(G + H)x 1 Hx 1 
—L0lX2 lx. 2 —1X2 


The characteristic polynomial is given by 


A 3 + /2A 2 + /iAT/o =0 (75) 

where 

h = Ix 2 + (G + H ) X1 + (E + uF)x 0 (76) 

fi = Axpx\ + EIX 0 X 2 + GIx 1 x 2 (77) 

fo = EGIxox\X 2 (78) 

and A = E(G + H) - (1 - uj)FG. 


To satisfy the Routh-Hurwitz criteria, we require fi > 0, 
for i = 0,1, 2 and / 1/2 — fo > 0. Notice that 

/ 1/2 - fo = (EIx 0 x 2 + Ax oXl )f 2 

+ (Ix 2 + (G + H)xi + ujFxo)GIxiX 2 (79) 

We always ha ve f n > 0 and /2 > 0. 

From Eq. |74| we deduce that A = (1 — oj)EG^j^ . There¬ 
fore, if F 3 (D) > 0, that is to say > 0, then A > 0. Hence, 
/1 > 0 and / 1/2 — fo > 0, so that SS3 is stable 

On the other hand, if < 0 and X 2 is very small, which 

occurs when SS3 is very close to SS2 11 , then /2 has the sign of 
A since the term with X 2 is negligible compared to the term 
Axoxi: 

f 2 = Ax 0 x 1 + (EIx 0 + xGIxi)x 2 < 0 
Thus, SS3 is unstable. □ 

Proof [Proposition[3] Since we always have fo > 0 and > 0, 
from the previous proof it follows that SS3 is stable if, and 
only if, / 1/2 — fo > 0. Indeed, this condition implies that we 
have also /1 > / 0//2 > 0. Using Eq. |49| and Eq. |79| we see 
that 

hf 2 - fo = F 4 (D,s i 0 n ) 

Therefore, the condition / 1/2 — fo > 0 is equivalent to 
F 4 (D, 4 n ) >0. □ 


D.3 Operating diagrams 

Proof [Proposition^] We know that SSI always exist and is 
stable. We know that SS2 b is unstable if it exists. Using Table 
[T] and Remark [ 2 ] we obtain the following results 

- fj\ is defined by D > D\ or 0 < D < D\ and S' c h,i n < 
F 4 (D)/YpY 4 . Therefore, SSI is the only existing steady 
state in this region. 

- J 2 if defined by D 3 < D < D± and Fi (0)/YoY 4 < S ch ,i n < 
F 2 (D)/Y 3 Y 4 . Therefore, both steady state SS2 exist and 
SS2# is stable since F 3 (D) > 0. 

- Jo if defined by 0 < D < D 2 and F 2 (D)/Y 3 Y 4 < 5 c h,in 
and F 4 (D, Sch.in/Ualii) > 0 when 0 < D < D 3 . Therefore, 
SS3 exists and is stable, both steady state SS2 exist and 
SS2# is unstable since F 3 (D) < 0.. 

- J 4 if defined by 0 < D < F 3 and Fi(D)/Y 3 Y 4 < f> c h,in < 
F 2 (D)/Y 3 Y 4 . Therefore, both steady state SS2 exist and 
SS2“ is unstable since F 3 (D) < 0. 

- Jh if defined by 0 < D < D 3 , F 2 (D)/Y 3 Y 4 < S ch ,i n and 
F 4 (D,S ctl , iD /Y 3 Y 4 ) < 0. Therefore, SS3 exists and is un¬ 
stable and both steady state SS2 exist and SS2** is unsta¬ 
ble since F 3 (D) < 0.. 

□ 

Proof [Propositions [ 5 ] [6] and [ 7 ] The result follows from Table 
^and Remark[2] The details are as in the proof of Proposition 

0 □ 


D.4 General case 

Proof [Lemma |4] Let S 2 > 0 be fixed. By H4, the function 

so e [0, + 00 ) >-> po(so, S 2 ) 6 [0, Po(Too, s 2 )) 

is monotonically increasing. Hence, it has an inverse function 
denoted by 


y S [0,po(+oo,S2)) s-t Mo(y,s 2 ) e [0,+oc), 

















Generalised approach to modelling a three-tiered microbial food-web 


17 


such that for all so > 0, s 2 > 0 and y £ [0, /to(+oo, s 2 )) ( |55[ > 
holds. □ 

Proof [LemmaJS] Let S 2 > 0 be fixed. By H5, the function 

si e [ 0 , + 00 ) >-> (si, S 2 ) e [0, /11 (+ 00 , S 2 )) 

is monotonically increasing. Hence, it has an inverse function 
denoted by 

y £ [0, /11 (+ 00 , S 2 )) m- Mi(y, S 2 ) S [0, + 00 ), 

such that for all si > 0, S 2 > 0 and y £ [0, £ [0, /11 (+ 00 , S 2 )) 
( |56| ) holds. □ 

Proof [Lemma lb] By H6, the function S 2 £ [0, + 00 ) m- 
y, 2 (S 2 ) S [0,M2(+oo)) is monotonically increasing. Hence, it 
has an inverse function denoted by 

y e [0,M2(+°o)) m- M 2 {y) £ [0,+oo), 

such that, for all S 2 > 0 and y £ [0,/t 2 (+oo)) ( |57| ) holds. □ 

Proof [LemmalT] By H7, for D-\-ao < /to(+oo, + 00 ) and D + 
a i < Mi(+oo,0jthere exist unique values s§ and Sj such that 
Eq. |58| holds, see Fig. [l|[a). □ 
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